What do we have in the training set?
csv = read.csv("../input/train.csv", header=T)
head(csv)
The description of the data says that “a number of Freesound audio samples were automatically annotated with labels from the AudioSet Ontology … Then, a data validation process was carried out in which a number of participants did listen to the annotated sounds and manually assessed the presence/absence of an automatically assigned sound category”. And that “The non-verified annotations of the train set have a quality estimate of at least 65-70% in each category.”
With that in mind, I would prefer to use only the manually verified data.
csv = csv[csv$manually_verified == 1,]
head(csv)
With that reduction, how many do we have in each category?
labalTable = table(csv$label)
There are 67 sound clips labeled “Bark”.
Of course, all of these files are in the inputs folder.
csv$filePath = paste0("../input/audio_train/", csv$fname)
Grab all of the sounds of dogs barking.
barkSoundClips = csv$filePath[csv$label=="Bark"]
library(audio)
## Warning: package 'audio' was built under R version 3.4.4
b1 = load.wave(barkSoundClips[1])
sample.rate = attr(b1, "rate")
str(b1)
## Class 'audioSample' atomic [1:620928] -3.05e-05 3.05e-05 0.00 -3.05e-05 -3.05e-05 ...
## ..- attr(*, "rate")= int 44100
## ..- attr(*, "bits")= int 16
The “rate” of 44,100 means that the pressure in front of the microphone was measured 44,100 times per second. Thus, the entire clip is 14 seconds long, and has 620,928 values.
play(b1)
<audio controls>
<source src="../input/audio_train/00e2b4cd.wav">
</audio>
plot(b1, type="l", col="darkblue", las=1)
start=45000
end=67000
abline(v=start, col="red"); abline(v=end, col="red")
Zoom in to the range marked by the two red lines.
plot(start:end, b1[start:end], type="l", col="darkblue", las=1)
abline(v=start, col="red"); abline(v=end, col="red")
play(b1[start:end])
arf = as.numeric(b1[start:end])
Lets pretend we sample some really ideal sound at a rate of 10 measurments per second. Let x be your time line in seconds.
example.rate = 10
x <- seq(-30, 70, by = 1/example.rate)
wav.freq = c(0.8, 1.5, 2)
wav.1 = sin(wav.freq[1] * 2 * pi * x)
wav.2 = sin(wav.freq[2] * 2 * pi * x)
wav.3 = rep(0, length(x))
xmid = seq(0, 20, by = 1/example.rate)
someMidXs = which(x>0 & x < max(x)/3)
xmid = x[someMidXs]
wav.3[someMidXs] = 1.6 * sin(wav.freq[3] * 2 * pi * xmid)
wav.sum = wav.1 + wav.2 + wav.3
waves = list("wave1" = wav.1,
"wave2" = wav.2,
"wave3" = wav.3,
"sum" = wav.sum)
The Forier transfrom is beautifully explained by 3blue1brown in his video “But what is the Fourier Transform? A visual introduction.” I’m using that as my bases for making the following function.
fourtrans = function(winding, g, t){
# winding - scaler, the "winding frequency"
# g - a vector describing the original wave.
# t - a vector of times corresponding to g; g and t must be the same length
i = complex(imaginary=1)
mass = g * exp( -2 * pi * i * winding * t)
return(sum(sum(mass))) # the double sum is diliberate.
}
Do the fourier transform on the example waves.
ex.freqs = seq(0.001, 3, 0.001)
sum.ft = sapply(ex.freqs, fourtrans, g=wav.sum, t=x)
plot(ex.freqs, abs(sum.ft), type="l", las=1, xlab="frequency")
polygon(x=c(max(ex.freqs), min(ex.freqs), ex.freqs), c(0,0, abs(sum.ft)), col="gray")
for (i in 1:length(wav.freq)){
abline(v=wav.freq[i], col=wav.colors[i])
text(x=wav.freq[i], y=par("usr")[4], xpd=T, pos=3, col=wav.colors[i], labels=wav.freq[i])
}
I see what I want. I see peaks that corresond to my original frequencies. But how to extract them when I don’t know the values and I don’t know how many there are? If I take every peak, I will have too many.
getMaxima = function(vector){
n=length(vector)
left = vector[1:(n-2)]
main = vector[2:(n-1)]
right = vector [3:n]
indecies = which(main > left & main > right)
# these indeces are for 'main' which is shifted from the original by 1.
return(indecies + 1)
}
rawMaxima = getMaxima(abs(sum.ft))
Using the raw values and taking EVERY maximum gives me 227 frequencies.
plot(ex.freqs, abs(sum.ft), type="l", las=1, xlab="frequency")
polygon(x=c(max(ex.freqs), min(ex.freqs), ex.freqs), c(0,0, abs(sum.ft)), col="gray")
points(col="red", pch=16, x=ex.freqs[rawMaxima], abs(sum.ft)[rawMaxima])
legend(legend=length(rawMaxima), x="topright", text.col="red", pch=16, col="red")
I could order them by their intensity (y-axis), but without knowing how many there are, I might stop with only 2 or take on several too many.
I know I can count on having a lot of little hills that don’t matter, so I’m going take the 75th percentile and use that a bench mark for “surely small”. Any maxima that are no higher than that, can be considered noise.
small = quantile(abs(sum.ft)[rawMaxima], .75)
dotColor = rep("red", length(rawMaxima))
dotColor[abs(sum.ft)[rawMaxima] <= small] = "lightblue"
plot(ex.freqs, abs(sum.ft), type="l", las=1, xlab="frequency")
polygon(x=c(max(ex.freqs), min(ex.freqs), ex.freqs), c(0,0, abs(sum.ft)), col="gray")
points(col=dotColor, pch=16, x=ex.freqs[rawMaxima], abs(sum.ft)[rawMaxima])
abline(h=small, lty=2, col="lightblue")
legend(legend=sum(dotColor=="red"), x="topright", text.col="red", pch=16, col="red")
Smooth out the values to remove unimportant peaks.
n.panels = 8
par(mfrow=c(n.panels,1), las=1, mar=c(0,3,0,0), oma=c(4,0,3,0))
plotSmoothMaxes = function(x, y, dotsAt, thresh){
#x - the x values of points to plot
#y - the y values of points to plot
#dotsAt - the indices of points to highlight (the maxes)
#thresh - the threshold y-value below which to ignore the maxes
plot(x, y, type="l", las=1, xlab="frequency", xaxt="n")
polygon(x=c(max(x), min(x), x), c(0,0, y), col="gray")
dotCol = rep("red", length(dotsAt))
dotCol[y[dotsAt] <= thresh] = "lightblue"
abline(h=thresh, lty=2, col="lightblue")
points(col=dotCol, pch=16, x=x[dotsAt], y[dotsAt])
}
span = seq(0.01, .2, length.out = 40)
n.span = length(span)
drawThese = round(seq(1, n.span, length.out = n.panels))
drawList = list()
numMaxes = rep(NA, length(span))
for (i in 1:n.span){
ft.lo = loess(abs(sum.ft) ~ ex.freqs, span=span[i])
smoother = predict(ft.lo, ex.freqs)
maxes = getMaxima(smoother)
numMaxes[i] = sum(smoother[maxes] > small)
if (i %in% drawThese) {
plotSmoothMaxes(x=ex.freqs, y=smoother, dotsAt = maxes, thresh=small)
abline(v=wav.freq, col=wav.colors)
legend(x="topright", legend=paste("span =", span[i]), bty="n")
}
}
title("smoothed data with different 'span'", outer = T)
Using a span of 0.07 is just enough to only get the main peaks. We still get the right number of peaks even with a span of more than twice that. As we increase the span, we see the values we would get from our maxima shift away from the accurate values (the red dots move away from the red, green and gold lines).
plot(span, numMaxes, type="l", las=1)
abline(h=3, lty=2, col="red")
Now we need a rule a computer can follow that will help it come to similar conclusions about optimizing the span for a given clip. There is probably some really good algorithm for doing that, but for now this will do:
# What is the most common number of peaks (across different span parameters)
# That must be right. What is the lowest span that gives me that many? That must be best.
getBestSpan = function(numMaxes, span, showplot=F){
tb = table(numMaxes)
tb = tb[order(tb, decreasing = T)]
mode = as.numeric(names(tb)[1])
w = which(numMaxes == mode)
lowSpan = span[min(w)]
if (showplot){
plot(span, numMaxes, type = "l", las=1)
abline(v=lowSpan, col="red")
abline(h=mode, col="blue", lty=2)
legend(x="topright", legend=paste("span =", lowSpan), bty="n")
title("Optimal Span")
}
return(lowSpan)
}
lowSpan = getBestSpan(numMaxes, span)
Based on this logic, 0.068 is the optimum span parameter for this clip.
getBestFreqs = function(ft, freqs, lowSpan, small, showplot=F){
ft.lo.best = loess(abs(ft) ~ freqs, span=lowSpan)
smoother = predict(ft.lo.best, freqs)
maxes = getMaxima(smoother)
maxes = maxes[smoother[maxes] > small]
freq.of.interest = freqs[maxes]
if(showplot){
plotSmoothMaxes(x=freqs, y=smoother, dotsAt = maxes, thresh=small)
}
return(freq.of.interest)
}
freq.of.interest = getBestFreqs(sum.ft, ex.freqs, lowSpan, small=small)
And based on that value for span, the frequencies of interest are:
## [1] 0.798 1.500 1.997
The peak for the wave 3 frequency is shorter and wider than the other two. This makes sense since it was only present in a small portion of the clip while the other two were present consitently. Lets apply a sliding window approach.
Define the windows.
winSize = 100 # how may values wide each window is
shift = 10 # how many values from the start of one window to the start of the next
len = length(wav.sum) # how many total values in the clip
winStarts = seq(1, len-winSize, shift)
winEnds = seq(winSize, len, shift)
nw = length(winStarts) # number of windows
For each window, do the fourier transform and save the value for each frequency.
#ex.freqs = seq(0.001, 1, 0.001)
shiftingFT = matrix(NA, ncol = length(ex.freqs), nrow=nw,
dimnames=list(startTime=x[winStarts], frequency=ex.freqs))
for (i in 1:nw){
win = wav.sum[winStarts[i]:winEnds[i]]
shiftingFT[i,] = sapply(ex.freqs, fourtrans, g=win, t=x[winStarts[i]:winEnds[i]])
}
shiftingFT.real = abs(shiftingFT)
midTime=round((x[winStarts]+x[winEnds])/2)
image(x=midTime, y=ex.freqs, z=shiftingFT.real, las=1, col=heat.colors(length(shiftingFT.real)))
For just the frequences of intersest, plot their intensity over time (across the sliding windows).
freq.of.interest = freq.of.interest [order(freq.of.interest)] # make sure they are in the original order
freqName = as.character(freq.of.interest)
plot(1,1, xlim=range(x), ylim=c(0,300), type="n",
xlab="time", ylab="frequency intensity", las=1)
for (i in 1:3){
fq = freqName[i]
lines(x=midTime, y=shiftingFT.real[,fq], col=wav.colors[i])
}
legend(x="topright", legend = round(freq.of.interest,3), col=wav.colors, lty=1)
The only thing we did to optimize the window size and shift size was… well nothing. I looked at the picture and adjusted it. I don’t know any good rules there, except hope that the same parameters work well over most cases.
Given this information, can you get the original data back? Because my measurements are for time intervales, and each interval must be >0, I don’t have values for the begining and end. We loose the edges of the data.
rebuilt = matrix(ncol=length(freq.of.interest),
nrow=length(x),
dimnames=list(x=x, freq=freqName))
for (i in 1:length(freq.of.interest)){
fq = freqName[i]
#.1 was the lowest span that I tried that didn't produce an error
intensityPerWindow.lo = loess(shiftingFT.real[,fq] ~ midTime, span=0.1)
amplitudeFactor = predict(intensityPerWindow.lo, x)
rebuilt[,i] = sin(freq.of.interest[i] * 2 * pi * x) * amplitudeFactor
}
Normalize the amplitude to roughly match the original.
maxWavSum = max(wav.sum)
intensitySumPerSec = rowSums(shiftingFT.real[,freqName])
maxIntesitySum = max(intensitySumPerSec)
ampFactor = maxWavSum / maxIntesitySum
rebuilt = rebuilt * ampFactor
Not perfect… but not too bad. I think real sound data will have more cycles per unit, and I think that will mean an increase in resolution. It would also look cleaner to have used a cutoff rather than scaling the amplitude by the intensity, but I don’t think that is as true for real data as it is for this simulation.
Here’s the original for comparison:
The original wave (g) is represented by our input data. For t, we need to create a vector of time measurements corresponding to the values in our sound clip.
getTimes = function(clipVals, rate=44100){
#clipVals - a vector representing a sound clip.
#rate - sample rate of the sound clip (probably 44,100)
# value => a vector representing the time measurements for those values
indecies = 1:length(clipVals)
return(indecies/rate)
}
times = getTimes(arf, rate=sample.rate)
Based on Wikipedia’s Audio Frequency page, sound frequencies that we can hear will be covered by a range of about 20 to 20,000 Hz.
freq = seq(20, 20000, 40)
Do the Fourier Transform for the arf sample.
arf.ft = sapply(freq, fourtrans, g=arf, t=times)
arfMaxes = getMaxima(abs(arf.ft))
# plot(freq, abs(arf.ft), type="l", las=1, xlab="frequency")
# title("Fourier transformed dog bark")
# polygon(x=c(max(freq), min(freq), freq), c(0,0, abs(arf.ft)), col="gray")
# points(col="red", pch=16, x=freq[arfMaxes], y=abs(arf.ft)[arfMaxes])
plotSmoothMaxes(x=freq, y=abs(arf.ft),
dotsAt = arfMaxes,
thresh=quantile(abs(arf.ft)[getMaxima(abs(arf.ft))], .95))
Which frequencies are of interest?
# Earlier I was interested in a lot of plotting.
# Here, just put the important bits into a function, no plots
getFreqOfInterest.old = function(ft, freqs,
small = quantile(abs(ft)[getMaxima(abs(ft))], .95),
span = seq(0.01, .2, length.out = 40)){
# ft = fourier transformed data
# ex.freqs - frequences to iterate over
# small = frequencies with intensities lower than this are discarded
# span - a vector of spans to iterate through
n.span = length(span)
numMaxes = rep(NA, length(span))
for (i in 1:n.span){
tryCatch(expr = {
numMaxes[i] = suppressWarnings(expr={
length(getBestFreqs(ft, freqs, lowSpan=span[i], small=small, showplot = T))
})
}, finally = {
#rm(ft.lo); rm(smoother); rm(maxes)
})
}
lowSpan = getBestSpan(numMaxes, span, showplot = T)
freq.of.interest = getBestFreqs(ft, freqs, lowSpan, small=small, showplot = T)
return(freq.of.interest)
}
getPeakFreqs <- function(ft, small=NULL){
if (is.null(small)){
fracMax = max(ft)/10
quant = quantile(abs(ft)[getMaxima(abs(ft))], .95)
small = max(fracMax, quant)
}
return(findpeaks(abs(ft), minpeakheight = small,
nups = 2, ndowns = 2, minpeakdistance = 5))
}
freqPeaks = getPeakFreqs(ft=abs(arf.ft))
arfFOI = freq[freqPeaks[,2]]
#arfFOI = getFreqOfInterest(ft=arf.ft, freqs=freq)
arfFOI
## [1] 860 500
For each window, do the fourier transform and save the value for each frequency.
getShiftingFT = function(rawSound, sampleRate,
winSize = 5000, shift = 1000,
freqs=freq){
# rawSound - values from a sound file
# sampleRate - number of values per second in rawSound
# winSize - number of values per window
# shift - number of values to shift over to start the next window
# freq - freqencies to measure
len = length(rawSound) # how many total values in the clip
times = getTimes(rawSound, rate = sampleRate)
winStarts = seq(1, len-winSize, shift)
winEnds = seq(winSize, len, shift)
midTime = (times[winStarts] + times[winEnds])/2
nw = length(winStarts) # number of windows
shiftingFT = matrix(NA, ncol = length(freqs), nrow=nw,
dimnames=list(midTime=round(midTime,3), frequency=freqs))
for (i in 1:nw){
win = rawSound[winStarts[i]:winEnds[i]]
shiftingFT[i,] = sapply(freqs, fourtrans, g=win, t=times[winStarts[i]:winEnds[i]])
}
shiftingFT.real = abs(shiftingFT)
return(shiftingFT.real)
}
arfMatrix = getShiftingFT(rawSound = arf, sampleRate = sample.rate)
midTime = as.numeric(rownames(arfMatrix))
image(x=midTime, y=freq, z=arfMatrix, las=1, col=heat.colors(length(arfMatrix)))
Plotting all of the frequencies shows a lot of blank space, just show the freqencies up to 2x the highest freq of interest.
maxFreqShow = max(arfFOI) * 1.2
w = which(freq <= maxFreqShow)
plotMatrix = arfMatrix[,w]
image(x=midTime, y=as.numeric(colnames(plotMatrix)), z=plotMatrix, las=1, col=heat.colors(length(plotMatrix)), ylab="frequency (Hz)")
#freq.of.interest = freq.of.interest [order(freq.of.interest)] # make sure they are in the original order
freqName = as.character(arfFOI)
arf.col = wav.colors # they can just happen to be the same
plot(1,1, xlim=range(midTime), ylim=c(0,max(arfMatrix)), type="n",
xlab="time", ylab="frequency intensity", las=1)
for (i in 1:length(arfFOI)){
fn = freqName[i]
lines(x=midTime, y=arfMatrix[,fn], col=arf.col[i])
}
legend(x="topright", legend = round(arfFOI,3), col=arf.col, lty=1)
Hear each frequency alone.
for (i in 1:length(arfFOI)){
toneTime = c(times, times+max(times), times+(2*max(times)))
s = sin(arfFOI[i] * 2 * pi * toneTime)
tone = audioSample(s, rate=sample.rate, bits=16)
play(tone)
wait(2)
}
Given this information, can you get the original data back? Because my measurements are for time intervales, and each interval must be >0, I don’t have values for the begining and end. We loose the edges of the data.
rebuilt = matrix(ncol=length(arfFOI),
nrow=length(times),
dimnames=list(time=times, freq=freqName))
for (i in 1:length(arfFOI)){
fq = freqName[i]
intensityPerWindow.lo = loess(c(0,arfMatrix[,fq],0) ~ c(min(times),midTime,max(times)))
amplitudeFactor = predict(intensityPerWindow.lo, times)
rebuilt[,i] = sin(arfFOI[i] * 2 * pi * times) * amplitudeFactor
}
Normalize the amplitude to roughly match the original.
maxWavSum = max(arf)
intensitySumPerSec = rowSums(arfMatrix[,freqName])
maxIntesitySum = max(intensitySumPerSec)
ampFactor = maxWavSum / maxIntesitySum
rebuilt = rebuilt * ampFactor
arfRebuilt = rowSums(rebuilt)
artArf = audioSample(arfRebuilt, rate=sample.rate, bits=16)
plot(start:end, arfRebuilt, type="l", col="darkblue", las=1)
abline(v=start, col="red"); abline(v=end, col="red")
play(artArf)
play(b1[start:end])
plot(start:end, b1[start:end], type="l", col="darkblue", las=1)
abline(v=start, col="red"); abline(v=end, col="red")
library("seewave")
## Warning: package 'seewave' was built under R version 3.4.4
tm = timer(b1,threshold=5,msmooth=c(50,0))
length(tm$s)
## [1] 264
This found 264 peaks, and I don’t see an easy way to merge adjascent peaks.
Smothing by maxes.
numB1 = abs(as.numeric(b1))
names(numB1) = getTimes(numB1)
nextTime = min(getTimes(numB1)) + max(getTimes(numB1))
vals = numB1
n = 9
naturalFloor = rep(NA, n)
inflatedShape = as.list(naturalFloor)
plot(x=getTimes(numB1), y=numB1, type="n", col="darkblue", las=1, xlim=c(1.1, 1.3))
for (i in 1:n){
vals = vals[getMaxima(vals)]
#tv = table(vals)
#tv = tv[order(tv, decreasing = T)]
naturalFloor[i] = min(vals)
inflatedShape[[i]] = vals
lines(x=as.numeric(names(vals)), y=vals, type="l", col=i)
#hist(vals, main=i, 200)
}
legend(x="topright", lty=1, col=1:n, legend=1:n)
#plot(x=as.numeric(names(b1m1)), y=b1m1, type="l", col="darkblue", xlab="seconds", las=1)
nextTimeLab = as.character(nextTime)
inflatedShape = lapply(inflatedShape, function(x){x1=c(0, x, 0); names(x1)[1]="0"; names(x1)[length(x1)]=nextTimeLab; return(x1)})
Collect amplitude envelopes.
numB1 = abs(as.numeric(b1))
timesB1 = getTimes(numB1)
xlim=NULL
plot(timesB1, numB1, type="l", col="darkblue", las=1, xlim=xlim)
maxAmp = max(numB1)
minpeakheight = 0.1 * maxAmp
threshold = 0.01 * maxAmp
#threshold = naturalFloor[3]
allEnv = findpeaks(numB1, nups=5, ndowns=5,
minpeakheight = minpeakheight,
threshold = threshold,
minpeakdistance = sample.rate * .1)
allEnv = data.frame(allEnv)
names(allEnv) = c("height", "peak", "start", "end")
allEnv$peakTime = getTimes(numB1)[allEnv$peak]
points(allEnv$peakTime, allEnv[,"height"], col="red", pch=16)
abline(v=getTimes(numB1)[allEnv$start], col="green")
abline(v=getTimes(numB1)[allEnv$end], col="red")
xlim=c(1.10, 1.5)
plot(timesB1, numB1, type="l", col="darkblue", las=1, xlim=xlim)
points(allEnv$peakTime, allEnv[,"height"], col="red", pch=16)
abline(v=getTimes(numB1)[allEnv$start], col="green")
abline(v=getTimes(numB1)[allEnv$end], col="red")
I wasn’t happy with the results from timer, and I didn’t like the output format. I was able to get reasonable peaks with findpeaks, but I didn’t like the start and end that was assigned to each peak. I want to use the peaks from findpeaks, but define the start and stop of each peak based on when the sound returns to some threshold.
To do this, I want to look for continuous areas of low values.
A value is only true if the surrounding n values are also true.
noLonelyTrue = function(bool, n=1, left=n, right=n, naEdges=NA){
# bool - a vectore of boolean values
# n - how many values on either side must be true for an original value to still be true
# left, right - number of values to the left (preceeding) or right (after) that must be true
end = length(bool)
# matrices are faster than loops.
# any given row in this matrix can be read as "is the original value true, how about the one before it, and before that..."
# the offset by one is acheived by adding a value at the end.
toTheLeft = matrix(data=rep(c(bool,NA),left+1)[1:(end * (left+1))], nrow=length(bool), ncol=left+1, byrow = F)
# same as above, but the other direction
toTheRight = matrix(data=rep(c(rev(bool), NA), right+1)[1:(end * (right+1))], nrow=length(bool), ncol=right+1, byrow = F)
# true values are treated as 1 in rowSums
newBool = rowSums(toTheLeft)==left+1 & rev(rowSums(toTheRight))==right+1
# the first and last n values don't have enough neighbors, by default they are NA values.
if (!is.na(naEdges)){
newBool[1:left] = naEdges
newBool[(length(newBool)-right+1):length(newBool)] = naEdges
}
return(newBool)
}
bool = c(T,T,F,T,T,T,T,T,T,T,T,F,F,T,T,T,T)
noLonelyTrue(bool, n=1)
## [1] NA FALSE FALSE FALSE TRUE TRUE TRUE TRUE TRUE TRUE FALSE
## [12] FALSE FALSE FALSE TRUE TRUE NA
Now define start and stop of peask based on when there is a continuous low value.
plot(timesB1, numB1, type="l", col="darkblue", las=1, xlim=xlim)
points(allEnv$peakTime, allEnv[,"height"], col="red", pch=16)
abline(h=threshold, col="pink")
# to test for values in pauses, use the level 3 smooth-by-maxima values
# shape.lo = loess(inflatedShape[[3]] ~ names(inflatedShape[[3]]))
# shape3 = predict(shape.lo, getTimes(numB1))
# threshold = max(threshold, min(shape3) * 5)
# lines(getTimes(numB1), shape3, col="orange", lwd=3, lty=2)
#which values are below threshold
peak=allEnv$peak
nrow=length(numB1)
ncol=length(peak)
bt1 = numB1 < threshold
bt2 = noLonelyTrue(bt1, n=20, naEdges = F)
bt = matrix(data=bt2, byrow = F, nrow=nrow, ncol=ncol)
beforePeak = rep(1:length(numB1), length(peak)) < rep(peak, each=length(numB1))
bp = matrix(data=beforePeak, byrow=F, nrow=nrow, ncol=ncol)
isBefore = bt & bp
startInd = apply(isBefore, 2, function(x){max(which(x))})
abline(v=getTimes(numB1)[startInd], col="green")
isAfter = bt & !bp
endInd = apply(isAfter, 2, function(x){min(which(x))})
abline(v=getTimes(numB1)[endInd], col="red")
xlim=NULL
plot(timesB1, numB1, type="l", col="darkblue", las=1, xlim=xlim)
points(allEnv$peakTime, allEnv[,"height"], col="red", pch=16)
abline(h=threshold, col="pink")
abline(v=getTimes(numB1)[startInd], col="green")
abline(v=getTimes(numB1)[endInd], col="red")